#########################################################################   
#                                 INFO                                  #   
#########################################################################  

# PROJECT: Gender Stereotypes, Political Leadership, and Voting Behavior in Tunisia
# PURPOSE: Main analysis
# CREATED: September 2019 by Alexandra Blackman & Marlette Jackson
# INPUTS:  tun_bjka.xlsx; tun_yg.xlsx
# OUTPUTS: main figures (amces for both surveys and vignette coefficients)

#########################################################################   
#                                 SETUP                                 #   
######################################################################### 

rm(list = ls()) 
setwd("~/Dropbox/Voters and Female Politicians/Political Behavior")


#########################################################################   
#                        LOAD PACKAGES                                  #   
######################################################################### 

need <- c("doBy","foreign", "lubridate","extrafont", "car", "plyr",
          "reshape", "openxlsx", "tidyr", "dplyr", "ggplot2", "ggthemes",
          "RColorBrewer", "rms", "lmtest", "sandwich", "msm",
          "stargazer","texreg","cowplot","plm","margins","prediction",
          "gplots","ggpubr", "sjmisc","sjPlot","gtable","cjoint","cregg",
          "scales","multiwayvcov") # packages needed
have <- need %in% rownames(installed.packages()) 
if(any(!have)) install.packages(need[!have]) 
invisible(lapply(need, library, character.only=T)) 

#########################################################################   
#                 LOAD CLUSTERING, F-TEST FUNCTIONS                     #   
#########################################################################   

super.cluster.fun <- function(model, cluster){
  vcovCL<-cluster.vcov(model, cluster)
  coef<-coeftest(model, vcovCL)
  return(coef)
}

f.test <- function(model, cluster){
  vcovCL<-cluster.vcov(model, cluster)
  ftest <- waldtest(model, vcov = vcovCL, test = "F")
  return(ftest)
}

#########################################################################   
#                       LOAD DATA                                       #   
#########################################################################   

## LOAD BJKA HOUSEHOLD SURVEY (JULY 2017)
tun_bjka <- read.xlsx("tun_bjka.xlsx")

## LOAD YOUGOV ONLINE SURVEY (APRIL 2019)
tun_yg <- read.xlsx("tun_yg.xlsx")


#########################################################################   
#                     BJKA HOUSEHOLD SURVEY                             #   
#########################################################################   

#########################################################################   
#                     BJKA: VIGNETTE EXPERIMENT                         #   
#########################################################################   

m_all <- lm(support ~ tassign_women, data=tun_bjka)
summary(m_all)

m_patr <- lm(support ~ tassign_women, data=tun_bjka[tun_bjka$patriarchal==1,])
summary(m_patr)

m_egal <- lm(support ~ tassign_women, data=tun_bjka[tun_bjka$patriarchal==0,])
summary(m_egal)

# Extract all coefficients 
bjkaFrame <- data.frame(var = rownames(summary(m_all)$coef),
                      coef = summary(m_all)$coef[, 1],
                      se = summary(m_all)$coef[, 2],
                      modelName = "Overall")

patrFrame <- data.frame(var = rownames(summary(m_patr)$coef),
                        coef = summary(m_patr)$coef[, 1],
                        se = summary(m_patr)$coef[, 2],
                        modelName = "Patriarchal")

egalFrame <- data.frame(var = rownames(summary(m_egal)$coef),
                        coef = summary(m_egal)$coef[, 1],
                        se = summary(m_egal)$coef[, 2],
                        modelName = "Egalitarian")

# merge data
allModelFrame <- data.frame(rbind(bjkaFrame, patrFrame, egalFrame))

# drop Intercept estimate
allModelFrame$drop <- ifelse(allModelFrame$var %in% 
                               grep("(Intercept)", allModelFrame$var, value=TRUE), 1, 
                             ifelse(is.na(allModelFrame$var), NA, 0))
allModelFrame <- allModelFrame[allModelFrame$drop==0,]

allModelFrame$var2 <-NA
allModelFrame$var2[allModelFrame$var=="tassign_womensecurity"] <- "Miriam Security"
allModelFrame$var2[allModelFrame$var=="tassign_womentraditional"] <- "Miriam Traditional"
allModelFrame$var2[allModelFrame$var=="tassign_womenwomen"] <- "Miriam Women"


allModelFrame[nrow(allModelFrame) + 1,] = list(var=NA, coef=0, se=0, modelName="Egalitarian",
                                               drop=0, var2="Miriam Control")

allModelFrame[nrow(allModelFrame) + 1,] = list(var=NA, coef=0, se=0, modelName="Patriarchal",
                                               drop=0, var2="Miriam Control")

allModelFrame[nrow(allModelFrame) + 1,] = list(var=NA, coef=0, se=0, modelName="Overall",
                                               drop=0, var2="Miriam Control")

allModelFrame$var2 <- factor(allModelFrame$var2, levels = c("Miriam Control",
                                                            "Miriam Women",
                                                            "Miriam Traditional",
                                                            "Miriam Security"))

allModelFrame$modelName <- factor(allModelFrame$modelName, levels = c("Egalitarian",
                                                                      "Patriarchal",
                                                                      "Overall"))

# Specify the width of confidence intervals
interval2 <- -qnorm((1-0.95)/2)  # 95% multiplier

# Plot
vignette_bjka <- ggplot(allModelFrame, aes(x=var2, y=coef, colour=modelName)) +
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) + 
  #coord_cartesian(ylim=c(-.5,1)) +
  ylim(-.5,1.2) +
  coord_flip() +
  geom_point(aes(x = var2, y = coef, shape = modelName, colour=modelName),
             size = 3, lwd = 1/2, position = position_dodge(width = 1/2)) +
  geom_errorbar(aes(ymin = coef - se*interval2, ymax = coef + se*interval2), width=.15,
                lwd = .75, position = position_dodge(width = 1/2)) +
  scale_shape_manual(values=c(16, 15, 17),name="Type", guide = guide_legend(reverse = TRUE)) +
  scale_color_manual(values=c('#56B4E9','#E69F00', '#999999'),name="Type", guide = guide_legend(reverse = TRUE)) +
  theme_bw() +
  theme(axis.text.x=element_text(size=12),
        axis.title.x=element_text(size=12),
        axis.text.y=element_text(size=10),
        axis.title.y=element_text(size=12),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank()) +
  xlab("Vignette Treatment") + 
  ylab("Coefficient")

vignette_bjka

png(file = "~/Dropbox/Voters and Female Politicians/Political Behavior/Figures/vignette_bjka.png",
    width = 6, height = 6, units = 'in', res = 300)
vignette_bjka
dev.off()

#########################################################################   
#                     BJKA: CONJOINT EXPERIMENT                         #   
#########################################################################   

tun_bjka$never_resp <- ifelse(is.na(tun_bjka$choice_1)==T &
                                is.na(tun_bjka$choice_3)==T &
                                is.na(tun_bjka$choice_5)==T &
                                is.na(tun_bjka$choice_7)==T, 1,0)

tun_bjka$sometimes_resp <- ifelse(is.na(tun_bjka$choice_1)==T |
                                    is.na(tun_bjka$choice_3)==T |
                                    is.na(tun_bjka$choice_5)==T |
                                    is.na(tun_bjka$choice_7)==T, 1,0)

tun_bjka$cjoint_response <- ifelse(tun_bjka$sometimes_resp==1,"Sometimes","Always")

tun_bjka$cjoint_response <- ifelse(tun_bjka$never_resp==1,"Never",tun_bjka$cjoint_response)

table(tun_bjka$cjoint_response)

tun_bjka_long <- reshape(data = tun_bjka, 
                         idvar = "responseid", 
                         varying = list(Sex=c("sex_1","sex_2","sex_3","sex_4","sex_5","sex_6","sex_7","sex_8"),
                                        Party=c("party_1","party_2","party_3","party_4","party_5","party_6","party_7","party_8"),
                                        Job=c("job_1","job_2","job_3","job_4","job_5","job_6","job_7","job_8"),
                                        Town=c("town_1","town_2","town_3","town_4","town_5","town_6","town_7","town_8"),
                                        Education=c("educ_1","educ_2","educ_3","educ_4","educ_5","educ_6","educ_7","educ_8"),
                                        Age=c("age_1","age_2","age_3","age_4","age_5","age_6","age_7","age_8"),
                                        Policy=c("policy_1","policy_2","policy_3","policy_4","policy_5","policy_6","policy_7","policy_8"),
                                        Experience=c("exper_1","exper_2","exper_3","exper_4","exper_5","exper_6","exper_7","exper_8"),
                                        choice=c("choice_1","choice_2","choice_3","choice_4","choice_5","choice_6","choice_7","choice_8"),
                                        contest=c("contest_1","contest_2","contest_3","contest_4","contest_5","contest_6","contest_7","contest_8")), 
                         direction="long", 
                         v.names = c("Sex","Party","Job","Town","Education","Age","Policy","Experience","choice","contest"),
                         sep="_")

tun_bjka_long_na2 <- tun_bjka_long %>% filter(!is.na(choice)) %>% filter(!is.na(patriarchal)) %>%
  mutate(Sex = factor(Sex, levels=c("Male","Female"))) %>%
  mutate(Party = factor(Party, levels=c("Initiative Party","Afek Tounes","UPL","Democratic Current",
                                        "Current of Love","al-Irada Movement","Popular Front",
                                        "Machrouu Tounes","Ennahda","Nidaa Tounes"))) %>%
  mutate(Job = factor(Job, levels = c("Director", "Lawyer","Teacher","Employee in the private sector","Political Activist","Union Activist"))) %>%
  mutate(Town = factor(Town, levels = c("within 10 km of here","within 20 km of here","within 30 km of here","within 40 km of here"))) %>%
  mutate(Education = factor(Education, levels = c("Bachelors","Masters","Studied outside of Tunisia"))) %>%
  mutate(Age = factor(Age, levels=c("70","60","50","40","30"))) %>%
  mutate(Policy = factor(Policy, levels=c("Improve security","Improve women's rights","Fight corruption","Improve the economic situation"))) %>%
  mutate(Experience = factor(Experience, levels=c("Has previously served in the government","Has not previously served in the government")))

tun_bjka_long_always <- tun_bjka_long_na2[tun_bjka_long_na2$cjoint_response=="Always",]

# load cregg package
library(cregg)

# save formula
f1 <- choice ~ Sex + Policy + Age + Education + Job + Experience + Town + Party

## Overall AMCE ##
amce_bjka1 <- amce(tun_bjka_long_always, f1, id = ~responseid, level_order="descending")
plot(amce_bjka1)

## Main AMCE with Overall, Patriarchal, and Egalitarian ##
amce_bjka2 <- amce(tun_bjka_long_always, f1, id = ~responseid, level_order="descending")

# amce for patriarchal vs. egalitarian
am_patr <- cj(tun_bjka_long_always, f1, id = ~ responseid, estimate = "amce", by = ~ type, level_order="descending")

## Create combined plot
amce_bjka2$type <- "Overall"
plot_vars <- c("statistic","feature", "level","estimate", "std.error",
               "z","p","lower","upper","type")

# overall amce
amce_all <- amce_bjka2[plot_vars]

# subgroup amce
amce_sub <- am_patr[plot_vars]
amce_bjka <- rbind(amce_all,amce_sub)

amce_bjka$type <- factor(amce_bjka$type,levels=c("Overall","Egalitarian","Patriarchal"))

amce_bjka_plot <- ggplot(amce_bjka, aes(x=level, y=estimate, colour=feature))+
  geom_hline(yintercept = 0, colour = gray(1/2), size=.25, linetype='dashed') + 
  coord_flip() +
#  ylim(-.3,.3) +
  geom_point(aes(x = level, y = estimate, colour=feature),
             size = 1, lwd = 1/2, position = position_dodge(width = 1/2)) +
  geom_errorbar(aes(ymin = lower, ymax = upper), width=.1,
                lwd = .5, position = position_dodge(width = 1/2)) + 
  facet_wrap(~type, ncol = 3L) +
  scale_x_discrete(limits = c("Nidaa Tounes","Ennahda","Machrouu Tounes","Popular Front",
                              "al-Irada Movement","Current of Love","Democratic Current",
                              "UPL","Afek Tounes","Initiative Party",
                              "within 40 km of here","within 30 km of here",
                              "within 20 km of here","within 10 km of here",
                              "Has not previously served in the government",
                              "Has previously served in the government",
                              "Union Activist", "Political Activist",
                              "Employee in the private sector","Teacher","Lawyer",
                              "Director","Studied outside of Tunisia","Masters",
                              "Bachelors","30","40","50","60","70",
                              "Improve the economic situation","Fight corruption",
                              "Improve women's rights","Improve security",
                              "Female","Male")) +
  theme_bw() +
  theme(axis.text.x=element_text(size=7),
        axis.title.x=element_text(size=11),
        axis.text.y=element_text(size=8),
        axis.title.y=element_text(size=11),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        legend.position="bottom",
        legend.title = element_blank()) +
  xlab("Attribute Level") + 
  ylab("Change in Pr(Candidate Preferred)")

amce_bjka_plot

png(file = "~/Dropbox/Voters and Female Politicians/Political Behavior/Figures/conjoint_amce_bjka.png",
    width = 8, height = 6, units = 'in', res = 300)
amce_bjka_plot
dev.off()


## Marginal Means + Difference for candidate gender ##
mm_patr <- cj(tun_bjka_long_always, choice ~ Sex, id = ~responseid, estimate = "mm", by = ~type, level_order="descending")
mm_diff_patr <- cj(tun_bjka_long_always, choice ~ Sex, id = ~responseid, estimate = "mm_diff", 
                   by = ~type, level_order="descending")

plot_vars <- c("BY","statistic","feature", "level","estimate", "std.error",
               "z","p","lower","upper")

mm_1 <- mm_patr[plot_vars]
mm_2 <- mm_diff_patr[plot_vars]
mm_bjka <- rbind(mm_1,mm_2)

mm_bjka$BY <- factor(mm_bjka$BY,levels=c("Patriarchal","Egalitarian","Patriarchal - Egalitarian"))
vline.dat <- data.frame(BY=levels(mm_bjka$BY), vl=c(0.5,0.5,0.0)) 
mm_plot<-plot(mm_bjka,vline=NULL) + ggplot2::geom_vline(aes(xintercept=vl), data=vline.dat,
                                                        size=.25, linetype='dashed') + 
  ggplot2::facet_wrap(~BY, ncol = 3L, scales = "free_x")  +
  scale_x_continuous(breaks = scales::pretty_breaks(5))

mm_plot


## Marginal Means + Difference for all variables ##
mm_patr2 <- cj(tun_bjka_long_always, f1, id = ~responseid, estimate = "mm", by = ~type, level_order="descending")
mm_diff_patr2 <- cj(tun_bjka_long_always, f1, id = ~responseid, estimate = "mm_diff", 
                    by = ~type, level_order="descending")

plot_vars <- c("BY","statistic","feature", "level","estimate", "std.error",
               "z","p","lower","upper")

mm_12 <- mm_patr2[plot_vars]
mm_22 <- mm_diff_patr2[plot_vars]
mm_bjka_all <- rbind(mm_12,mm_22)

mm_bjka_all$BY <- factor(mm_bjka_all$BY,levels=c("Patriarchal","Egalitarian","Patriarchal - Egalitarian"))
vline.dat <- data.frame(BY=levels(mm_bjka_all$BY), vl=c(0.5,0.5,0.0)) 
mm_plot2<-plot(mm_bjka_all,vline=NULL) + ggplot2::geom_vline(aes(xintercept=vl), data=vline.dat,
                                                             size=.25, linetype='dashed') + 
  ggplot2::facet_wrap(~BY, ncol = 3L, scales = "free_x")  +
  scale_x_continuous(breaks = scales::pretty_breaks(5)) +
  ggplot2::theme(axis.text.x=element_text(size=7),
        axis.title.x=element_text(size=11),
        axis.text.y=element_text(size=8),
        axis.title.y=element_text(size=11),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        legend.position="bottom",
        legend.title = element_blank())

mm_plot2

png(file = "~/Dropbox/Voters and Female Politicians/Political Behavior/Figures/conjoint_mm_bjka.png",
    width = 8, height = 6, units = 'in', res = 300)
mm_plot2
dev.off()

#########################################################################   
#                   YOUGOV ONLINE SURVEY                                #   
#########################################################################   

#########################################################################   
#                   YOUGOV: VIGNETTE EXPERIMENT                         #   
#########################################################################   

tun_yg_long <- gather(tun_yg, condition, support, woman_control:man_security, factor_key=TRUE)

m_all <- lm(support ~ condition + respondentid, data=tun_yg_long)
summary(m_all)

m_patr <- lm(support ~ condition + respondentid, data=tun_yg_long[tun_yg_long$patriarchal==1,])
summary(m_patr)

m_egal <- lm(support ~ condition + respondentid, data=tun_yg_long[tun_yg_long$patriarchal==0,])
summary(m_egal)

# Cluster by respondent
overall_clust <- super.cluster.fun(m_all, tun_yg_long$respondentid)
patr_clust <- super.cluster.fun(m_patr, tun_yg_long$respondentid[tun_yg_long$patriarchal==1])
egal_clust <- super.cluster.fun(m_egal, tun_yg_long$respondentid[tun_yg_long$patriarchal==0])

# Extract coefficients + se from cluster object
ygFrame <- data.frame(var = rownames(summary(m_all)$coef),
                      coef = overall_clust[, 1],
                      se = overall_clust[, 2],
                      modelName = "Overall")

patrFrame <- data.frame(var = rownames(summary(m_patr)$coef),
                        coef = patr_clust[, 1],
                        se = patr_clust[, 2],
                        modelName = "Patriarchal")

egalFrame <- data.frame(var = rownames(summary(m_egal)$coef),
                        coef = egal_clust[, 1],
                        se = egal_clust[, 2],
                        modelName = "Egalitarian")

# merge data
allModelFrame <- data.frame(rbind(egalFrame, patrFrame, ygFrame))

# Drop intercept and fixed effect coefficients
allModelFrame$drop <- ifelse(allModelFrame$var %in% 
                               grep("respondentid", allModelFrame$var, value=TRUE), 1, 
                             ifelse(is.na(allModelFrame$var), NA, 0))

allModelFrame <- allModelFrame[allModelFrame$drop==0,]

allModelFrame$drop <- ifelse(allModelFrame$var %in% 
                               grep("(Intercept)", allModelFrame$var, value=TRUE), 1, 
                             ifelse(is.na(allModelFrame$var), NA, 0))

allModelFrame <- allModelFrame[allModelFrame$drop==0,]


allModelFrame$var2 <-NA
allModelFrame$var2[allModelFrame$var=="conditionman_control"] <- "Ahmed Control"
allModelFrame$var2[allModelFrame$var=="conditionwoman_security"] <- "Miriam Security"
allModelFrame$var2[allModelFrame$var=="conditionman_security"] <- "Ahmed Security"


allModelFrame[nrow(allModelFrame) + 1,] = list(var=NA, coef=0, se=0, modelName="Egalitarian",
                                               drop=0, var2="Miriam Control")

allModelFrame[nrow(allModelFrame) + 1,] = list(var=NA, coef=0, se=0, modelName="Patriarchal",
                                               drop=0, var2="Miriam Control")

allModelFrame[nrow(allModelFrame) + 1,] = list(var=NA, coef=0, se=0, modelName="Overall",
                                               drop=0, var2="Miriam Control")

allModelFrame$var2 <- factor(allModelFrame$var2, levels = c("Miriam Control",
                                                            "Ahmed Control",
                                                            "Miriam Security",
                                                            "Ahmed Security"))

allModelFrame$modelName <- factor(allModelFrame$modelName, levels = c("Egalitarian",
                                                                      "Patriarchal",
                                                                      "Overall"))


# Specify the width of confidence intervals
interval2 <- -qnorm((1-0.95)/2)  # 95% multiplier

# Plot

vignette_yg <- ggplot(allModelFrame, aes(x=var2, y=coef, colour=modelName)) +
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
  #coord_cartesian(ylim=c(-.5,1)) +
  ylim(-.5,1.2) +
  coord_flip() +
  geom_point(aes(x = var2, y = coef, shape = modelName, colour=modelName),
             size = 3, lwd = 1/2, position = position_dodge(width = 1/2)) +
  geom_errorbar(aes(ymin = coef - se*interval2, ymax = coef + se*interval2), width=.15,
                lwd = .75, position = position_dodge(width = 1/2)) +
  scale_shape_manual(values=c(16, 15, 17),name="Type", guide = guide_legend(reverse = TRUE))+
  scale_color_manual(values=c('#56B4E9','#E69F00', '#999999'),name="Type", guide = guide_legend(reverse = TRUE))+
  theme_bw() +
  theme(axis.text.x=element_text(size=12),
        axis.title.x=element_text(size=12),
        axis.text.y=element_text(size=10),
        axis.title.y=element_text(size=12),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank()) +
  #ggtitle("Vignette effects") + 
  xlab("Vignette Treatment") + 
  ylab("Coefficient")

vignette_yg

png(file = "~/Dropbox/Voters and Female Politicians/Political Behavior/Figures/vignette_yg.png",
    width = 6, height = 6, units = 'in', res = 300)
vignette_yg
dev.off()


#########################################################################   
#                   YOUGOV: CONJOINT EXPERIMENT                         #   
#########################################################################   

tun_yg$choice_profileA <- ifelse(tun_yg$QS3_cj1=="Candidate A",1,0)
tun_yg$choice_profileB <- ifelse(tun_yg$QS3_cj1=="Candidate B",1,0)

tun_yg$choice_profileC <- ifelse(tun_yg$QS3_cj2=="Candidate C",1,0)
tun_yg$choice_profileD <- ifelse(tun_yg$QS3_cj2=="Candidate D",1,0)

tun_yg$choice_profileE <- ifelse(tun_yg$QS3_cj3=="Candidate E",1,0)
tun_yg$choice_profileF <- ifelse(tun_yg$QS3_cj3=="Candidate F",1,0)

tun_yg$choice_profileG <- ifelse(tun_yg$QS3_cj4=="Candidate G",1,0)
tun_yg$choice_profileH <- ifelse(tun_yg$QS3_cj4=="Candidate H",1,0)

tun_yg$choice_profileI <- ifelse(tun_yg$QS3_cj5=="Candidate I",1,0)
tun_yg$choice_profileJ <- ifelse(tun_yg$QS3_cj5=="Candidate J",1,0)


tun_yg_conjoint <- reshape(data = tun_yg,
                           idvar = "RecordNo", 
                           varying = list(Sex=c("A1_profileA","A1_profileB","A1_profileC","A1_profileD","A1_profileE","A1_profileF",
                                                "A1_profileG","A1_profileH","A1_profileI","A1_profileJ"),
                                          Party=c("A2_profileA","A2_profileB","A2_profileC","A2_profileD","A2_profileE","A2_profileF",
                                                  "A2_profileG","A2_profileH","A2_profileI","A2_profileJ"),
                                          Job=c("A3_profileA","A3_profileB","A3_profileC","A3_profileD","A3_profileE","A3_profileF",
                                                "A3_profileG","A3_profileH","A3_profileI","A3_profileJ"),
                                          Town=c("A4_profileA","A4_profileB","A4_profileC","A4_profileD","A4_profileE","A4_profileF",
                                                 "A4_profileG","A4_profileH","A4_profileI","A4_profileJ"),
                                          Education=c("A5_profileA","A5_profileB","A5_profileC","A5_profileD","A5_profileE","A5_profileF",
                                                      "A5_profileG","A5_profileH","A5_profileI","A5_profileJ"),
                                          Experience=c("A6_profileA","A6_profileB","A6_profileC","A6_profileD","A6_profileE","A6_profileF",
                                                       "A6_profileG","A6_profileH","A6_profileI","A6_profileJ"),
                                          Policy=c("A7_profileA","A7_profileB","A7_profileC","A7_profileD","A7_profileE","A7_profileF",
                                                   "A7_profileG","A7_profileH","A7_profileI","A7_profileJ"),
                                          Age=c("A8_profileA","A8_profileB","A8_profileC","A8_profileD","A8_profileE","A8_profileF",
                                                "A8_profileG","A8_profileH","A8_profileI","A8_profileJ"),
                                          choice=c("choice_profileA","choice_profileB","choice_profileC","choice_profileD","choice_profileE",
                                                   "choice_profileF","choice_profileG","choice_profileH","choice_profileI","choice_profileJ")), 
                           direction="long", 
                           v.names = c("Sex","Party","Job","Town","Education","Experience","Policy","Age","choice"),
                           sep="_")

table(tun_yg_conjoint$Party)

# Rename Jabha Shabiyya --> Popular Front
tun_yg_conjoint$Party <- ifelse(tun_yg_conjoint$Party=="Jabha Shabiyya", "Popular Front",
                                 tun_yg_conjoint$Party)

# Rename Studied outside of Tunis --> Studied outside of Tunisia
tun_yg_conjoint$Education <- ifelse(tun_yg_conjoint$Education=="Studied outside of Tunis", "Studied outside of Tunisia",
                                tun_yg_conjoint$Education)

# Rename Director of an organization --> Director
tun_yg_conjoint$Job <- ifelse(tun_yg_conjoint$Job=="Director of an organization", "Director",
                                    tun_yg_conjoint$Job)


tun_yg_conjoint_na <- tun_yg_conjoint %>% filter(!is.na(choice)) %>% filter(!is.na(patriarchal)) %>%
  mutate(Sex = factor(Sex, levels=c("Male","Female"))) %>%
  mutate(Party = factor(Party, levels=c("Independent List","Democratic Current","Popular Front",
                                        "Ennahda","Nidaa Tounes"))) %>%
  mutate(Job = factor(Job, levels = c("Director", "Lawyer","Teacher",
                                      "Private sector employee","Political Activist","Union Activist"))) %>%
  mutate(Town = factor(Town, levels = c("Within 10 km of your municipality",
                                        "Within 20 km of your municipality",
                                        "Within 30 km of your municipality",
                                        "Within 40 km of your municipality"))) %>%
  mutate(Education = factor(Education, levels = c("Bachelors","Masters","Studied outside of Tunisia"))) %>%
  mutate(Age = factor(Age, levels=c("71","60","51","42","33"))) %>%
  mutate(Policy = factor(Policy, levels=c("Improve security","Improve women's rights","Fight corruption","Improve the economic situation"))) %>%
  mutate(Experience = factor(Experience, levels=c("Has previously served in the government","Has not previously served in the government")))



library(cregg)

# save formula
f1 <- choice ~ Sex + Policy + Age + Education + Job + Experience + Town + Party

## Overall AMCE ##
amce_yg1 <- amce(tun_yg_conjoint_na, f1, id = ~RecordNo, level_order="descending")
plot(amce_yg1)

## Main AMCE with Overall, Patriarchal, and Egalitarian ##
# overall amce
amce_yg2 <- amce(tun_yg_conjoint_na, f1, id = ~RecordNo, level_order="descending")

# amce for patriarchal vs. egalitarian
am_patr <- cj(tun_yg_conjoint_na, f1, id = ~ RecordNo, estimate = "amce", by = ~ type, level_order="descending")

# plot all three together
amce_yg2$type <- "Overall"
plot_vars <- c("statistic","feature", "level","estimate", "std.error",
               "z","p","lower","upper","type")

amce_all <- amce_yg2[plot_vars]
amce_sub <- am_patr[plot_vars]
amce_yg <- rbind(amce_all,amce_sub)


amce_yg$type <- factor(amce_yg$type,levels=c("Overall","Egalitarian","Patriarchal"))

amce_yg_plot <- ggplot(amce_yg, aes(x=level, y=estimate, colour=feature))+
  geom_hline(yintercept = 0, colour = gray(1/2), size=.25, linetype='dashed') + 
  coord_flip() +
  #  ylim(-.3,.3) +
  geom_point(aes(x = level, y = estimate, colour=feature),
             size = 1, lwd = 1/2, position = position_dodge(width = 1/2)) +
  geom_errorbar(aes(ymin = lower, ymax = upper), width=.1,
                lwd = .5, position = position_dodge(width = 1/2)) + 
  facet_wrap(~type, ncol = 3L) +
  scale_x_discrete(limits = c("Nidaa Tounes","Ennahda","Popular Front",
                              "Democratic Current","Independent List",
                              "Within 40 km of your municipality",
                              "Within 30 km of your municipality",
                              "Within 20 km of your municipality",
                              "Within 10 km of your municipality",
                              "Has not previously served in the government",
                              "Has previously served in the government",
                              "Union Activist", "Political Activist",
                              "Private sector employee","Teacher","Lawyer",
                              "Director","Studied outside of Tunisia",
                              "Masters","Bachelors","33","42","51","60","71",
                              "Improve the economic situation","Fight corruption",
                              "Improve women's rights","Improve security",
                              "Female","Male")) +
  theme_bw() +
  theme(axis.text.x=element_text(size=7),
        axis.title.x=element_text(size=11),
        axis.text.y=element_text(size=8),
        axis.title.y=element_text(size=11),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        legend.position="bottom",
        legend.title = element_blank()) +
  xlab("Attribute Level") + 
  ylab("Change in Pr(Candidate Preferred)")

amce_yg_plot

png(file = "~/Dropbox/Voters and Female Politicians/Political Behavior/Figures/conjoint_amce_yg.png",
    width = 8, height = 6, units = 'in', res = 300)
amce_yg_plot
dev.off()


## Marginal Means + Difference for candidate gender ##
mm_patr <- cj(tun_yg_conjoint_na, choice ~ Sex, id = ~RecordNo, estimate = "mm", by = ~type, level_order="descending")
mm_diff_patr <- cj(tun_yg_conjoint_na, choice ~ Sex, id = ~RecordNo, estimate = "mm_diff", 
                   by = ~type, level_order="descending")

plot_vars <- c("BY","statistic","feature", "level","estimate", "std.error",
               "z","p","lower","upper")

mm_1 <- mm_patr[plot_vars]
mm_2 <- mm_diff_patr[plot_vars]
mm_yg <- rbind(mm_1,mm_2)

mm_yg$BY <- factor(mm_yg$BY,levels=c("Patriarchal","Egalitarian","Patriarchal - Egalitarian"))
vline.dat <- data.frame(BY=levels(mm_yg$BY), vl=c(0.5,0.5,0.0)) 

my_breaks <- function(x) { if (max(x) < .2) seq(-.2, .2, .2) else seq(.4, .6, .1) }


mm_plot<-plot(mm_yg,vline=NULL) + ggplot2::geom_vline(aes(xintercept=vl), data=vline.dat,
                                                      size=.25, linetype='dashed') + 
  ggplot2::facet_wrap(~BY, ncol = 3L, scales = "free_x")  +
  scale_x_continuous(breaks = scales::pretty_breaks(5))

mm_plot


## Marginal Means + Difference for all variables ##
mm_patr2 <- cj(tun_yg_conjoint_na, f1, id = ~RecordNo, estimate = "mm", by = ~type, level_order="descending")
mm_diff_patr2 <- cj(tun_yg_conjoint_na, f1, id = ~RecordNo, estimate = "mm_diff", 
                    by = ~type, level_order="descending")

plot_vars <- c("BY","statistic","feature", "level","estimate", "std.error",
               "z","p","lower","upper")

mm_12 <- mm_patr2[plot_vars]
mm_22 <- mm_diff_patr2[plot_vars]
mm_yg_all <- rbind(mm_12,mm_22)

mm_yg_all$BY <- factor(mm_yg_all$BY,levels=c("Patriarchal","Egalitarian","Patriarchal - Egalitarian"))
vline.dat <- data.frame(BY=levels(mm_yg_all$BY), vl=c(0.5,0.5,0.0)) 

my_breaks <- function(q) { if (max(q) < .3) seq(-.1, .2, .1) else seq(.3, .7, .1) }

mm_plot2<-plot(mm_yg_all,vline=NULL) + ggplot2::geom_vline(aes(xintercept=vl), data=vline.dat,
                                                           size=.25, linetype='dashed') + 
  ggplot2::facet_wrap(~BY, ncol = 3L, scales = "free_x")  +
  #scale_x_continuous(breaks = scales::pretty_breaks(5))
  scale_x_continuous(breaks = my_breaks)+
  ggplot2::theme(axis.text.x=element_text(size=7),
                 axis.title.x=element_text(size=11),
                 axis.text.y=element_text(size=8),
                 axis.title.y=element_text(size=11),
                 panel.grid.major = element_blank(),
                 panel.grid.minor = element_blank(),
                 panel.background = element_blank(),
                 legend.position="bottom",
                 legend.title = element_blank())

mm_plot2

png(file = "~/Dropbox/Voters and Female Politicians/Political Behavior/Figures/conjoint_mm_yg.png",
    width = 8, height = 6, units = 'in', res = 300)
mm_plot2
dev.off()

